{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5c2f94e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd \n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.stats import norm\n",
    "from scipy.optimize import minimize\n",
    "from scipy.optimize import bisect"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a1ac26d",
   "metadata": {},
   "source": [
    "# Model\n",
    "\n",
    "The model is\n",
    "$$\n",
    "\\max_{m_i} \\frac{((1-\\tau)y_i-\\mu pm_i)^{1-\\sigma}}{1-\\sigma} + \\phi (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "The FOC is \n",
    "$$\n",
    "\\mu p ((1-\\tau)y_i- \\mu pm_i)^{-\\sigma} = \\phi\\alpha_0 \\exp(-\\alpha_1 - \\alpha_0 m_i)\n",
    "$$\n",
    "This gives a solution $m_i = {\\cal M}(y_i,p|\\sigma, \\phi, \\alpha_0, \\alpha_1, \\mu)$\n",
    "\n",
    "Health insurance system\n",
    "$$\n",
    "\\tau \\sum_i \\omega_i y_i = (1-\\mu)p \\sum_i \\omega_i m_i~~~\\Rightarrow~~~\\tau = (1-\\mu)s\n",
    "$$\n",
    "where $s=pm/y$ with $m=\\sum_i \\omega_i m_i$ and $y=\\sum_i \\omega_i y_i$\n",
    "\n",
    "\n",
    "### US\n",
    "For the US, we have $p=1$ and $\\sum_i \\omega_i y_i \\equiv y = 1$. Therefore, the moments restrictions  \n",
    "$$\n",
    "\\Psi_1 = \\frac{p \\sum_i \\omega_i m_i}{\\sum_i \\omega_i y_i} = \\sum_i \\omega_i m_i\n",
    "$$\n",
    "$$\n",
    "\\Psi_2 = \\sum_i \\omega_i (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "$$\n",
    "\\Psi_3 = \\frac{(1-\\exp(-\\alpha_1 - \\alpha_0 m_{\\max}))}{(1-\\exp(-\\alpha_1 - \\alpha_0 m_{\\min}))}\n",
    "$$\n",
    "allows to identify $\\{\\phi,\\alpha_0,\\alpha_1\\}$ with $\\tau = (1-\\mu) \\Psi_1$\n",
    "\n",
    "\n",
    "### Europe\n",
    "\n",
    "For Europe, we keep $\\{\\phi,\\alpha_0,\\alpha_1\\}$ found in the previous step, and the moments restrictions\n",
    "$$\n",
    "\\Psi_1 = \\frac{p \\sum_i \\omega_i m_i}{\\sum_i \\omega_i y_i}\n",
    "$$\n",
    "$$\n",
    "\\Psi_2 = \\sum_i \\omega_i (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "allows to identify $\\{p,\\alpha_1\\}$, given that a calibration of $\\sum_i \\omega_i y_i \\equiv y$ and $\\tau = (1-\\mu) \\Psi_1$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e15e6936",
   "metadata": {},
   "source": [
    "We start by programming functions for finding optimal $m$ from the first-order condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "01fd442f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def foc(m,sigma_c,phi,alphas,y,p,mu,mom):\n",
    "\treturn mu*p*( (1-(1-mu)*mom)*y - mu*p*m )**(-sigma_c) - alphas[0]*phi*np.exp(-alphas[1]-alphas[0]*m)\n",
    "\n",
    "def u(m,sigma_c,phi,alphas,y,p,mu,mom):\n",
    "\treturn ( ((1-(1-mu)*mom)*y - mu*p*m)**(1-sigma_c) )/(1-sigma_c) + phi*(1-np.exp(-alphas[1] - alphas[0]*m))\n",
    "\n",
    "def optm(sigma_c,phi,alphas,y,p,mu,mom):\n",
    "\ttry:\n",
    "\t\tm = bisect(foc,0.0,y*0.99,args=(sigma_c,phi,alphas,y,p,mu,mom))\n",
    "\texcept:\n",
    "\t\tup  = u(0.99*y,sigma_c,phi,alphas,y,p,mu,mom)\n",
    "\t\tlow = u(0.0,sigma_c,phi,alphas,y,p,mu,mom)\n",
    "\t\tif up > low:\n",
    "\t\t\tm = 0.99*y\n",
    "\t\telse :\n",
    "\t\t\tm = 0.0\n",
    "\treturn m\n",
    "\n",
    "def prod(alphas,m):\n",
    "\treturn 1 - np.exp(-alphas[1] - alphas[0]*m)\n",
    "\n",
    "def hopt(sigma_c,phi,alphas,y,p):\n",
    "\tm = optm(sigma_c,phi,alphas,y,p,mu,mom)\n",
    "\treturn prod(alphas,m)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87cca3d1",
   "metadata": {},
   "source": [
    "These are the parameters that we end up getting, but we check here that the function works. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "72da1031",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma_c   = 2\n",
    "p_us      = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1273ad26",
   "metadata": {},
   "source": [
    "We load up the parameters of the income process and use that to compute 4 median incomes within quartiles. The stationary variance of the income process of each country is "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8a9b58b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigs = [0.256367, 0.094643, 0.214538, 0.237439, 0.169834, 0.094643, 0.270357, 0.486429]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21276dcb",
   "metadata": {},
   "source": [
    "For Europe, we take a mean. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "54c54890",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.486429, 0.19111728571428574)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sig_us = sigs[-1]\n",
    "sig_eu = np.mean(sigs[0:7])\n",
    "sig_us,sig_eu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5096d3a0",
   "metadata": {},
   "source": [
    "We program up the income distribution, making sure it has a mean of 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "38b7dad4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.37922262, 0.67735473, 1.05644178, 1.88698087]), 1.0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qs    = [0.125,0.375,0.625,0.875]\n",
    "yy_us = [norm(0,np.sqrt(sig_us)).ppf(q) for q in qs]\n",
    "yy_us = np.exp(yy_us)\n",
    "yy_us = yy_us/np.mean(yy_us)\n",
    "yy_us,np.mean(yy_us)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9ce6034",
   "metadata": {},
   "source": [
    "We do the same for Europe but scale down everything to have an average income of 0.78"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5f4c9032",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.44109972, 0.63452033, 0.83837729, 1.20600267]), 0.78)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_eu  = 0.78\n",
    "qs    = [0.125,0.375,0.625,0.875]\n",
    "yy_eu = [norm(0,np.sqrt(sig_eu)).ppf(q) for q in qs]\n",
    "yy_eu = np.exp(yy_eu)\n",
    "yy_eu = yy_eu/np.mean(yy_eu)\n",
    "yy_eu = y_eu*yy_eu\n",
    "yy_eu,np.mean(yy_eu)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e279f29",
   "metadata": {},
   "source": [
    "We use the following moments for the U.S. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5bafd2b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "s_us    = 0.1504\n",
    "h_us    = 0.891\n",
    "grad_us = 1.276\n",
    "mu_us   = 0.136"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "388ef5ca",
   "metadata": {},
   "source": [
    "We use the following moments for Europe (we do not use the gradient but still compare):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b139b8ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "s_eu    = 0.0955\n",
    "h_eu    = 0.931\n",
    "grad_eu = 1.06\n",
    "mu_eu   = 0.155"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4ece99cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "no_social_security = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "14de71ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "if no_social_security:\n",
    "    mu_us = 1\n",
    "    mu_eu = 1\n",
    "else:\n",
    "    mu_us   = 0.136     \n",
    "    mu_eu   = 0.155"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e52353c8",
   "metadata": {},
   "source": [
    "# Calibration\n",
    "\n",
    "We start the exercise with finding U.S. parameters that fit the moments. We are solving a set of equations for three unknowns, $\\alpha_0$, $\\alpha_{1,US}$ and $\\phi$. We use the methods of moments to minimize the squared distance of moments. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a9f55003",
   "metadata": {},
   "outputs": [],
   "source": [
    "def criterion_us(theta,sigma_c,alphas,y,p,mu,moms):\n",
    "\talphas[0] = theta[0]\n",
    "\talphas[1] = theta[1]\n",
    "\tphi = theta[2]\n",
    "\tms = [optm(sigma_c,phi,alphas,i,p,mu,moms[0]) for i in y]\n",
    "\tmoms_s = np.zeros(3)\n",
    "\tmoms_s[0] = p*np.sum(ms)/np.sum(y)\n",
    "\tmoms_s[1] = np.mean([prod(alphas,m) for m in ms])\n",
    "\tmoms_s[2] = prod(alphas,ms[3])/prod(alphas,ms[0])\n",
    "\treturn (moms[0] - moms_s[0])**2 + (moms[1] - moms_s[1])**2 + (moms[2] - moms_s[2])**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b05de664",
   "metadata": {},
   "source": [
    "These are the estimates reported in the text:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f419d91b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([7.5500684 , 1.46959713]), 0.4043119502544598)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alphas_us    = np.zeros(2)\n",
    "opt_us       = minimize(criterion_us,[5.0,1.5,2],args=(sigma_c,[0.0,0.0],yy_us,p_us,mu_us,[s_us,h_us,grad_us]))\n",
    "alphas_us[0] = opt_us.x[0]\n",
    "alphas_us[1] = opt_us.x[1]\n",
    "phi          = opt_us.x[2]\n",
    "alphas_us,phi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2c30303",
   "metadata": {},
   "source": [
    "We can check that we match moments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "02e9c4af",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.15040831544413752, 0.8909971567754102, 1.275998560456296)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ms_us = [optm(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us) for i in yy_us]\n",
    "np.mean(ms_us), np.mean([prod(alphas_us,m) for m in ms_us]), prod(alphas_us,ms_us[3])/prod(alphas_us,ms_us[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "64e69769",
   "metadata": {},
   "outputs": [],
   "source": [
    "hs_us = [prod(alphas_us,m) for m in ms_us]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5eacecdf",
   "metadata": {},
   "source": [
    "In Europe, we seek to estimate $p_{EU}$ and $\\alpha_{1,EU}$ using $(s_{EU},h_{EU})$ as moments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "bcf770d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def criterion_eu(theta,sigma_c,phi,alphas,y,mu,moms):\n",
    "\talphas[1] = theta[0]\n",
    "\tp = theta[1]\n",
    "\tms = [optm(sigma_c,phi,alphas,i,p,mu,moms[0]) for i in y]\n",
    "\tmoms_s = np.zeros(2)\n",
    "\tmoms_s[0] = p*np.mean(ms)/np.mean(y)\n",
    "\tmoms_s[1] = np.mean([prod(alphas,m) for m in ms])\n",
    "\treturn (moms[0] - moms_s[0])**2 + (moms[1] - moms_s[1])**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d05ed9b",
   "metadata": {},
   "source": [
    "These are the estimates reported in the paper:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "4936d7b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.7067741235304366, 0.46200735305328194)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alphas_eu    = np.zeros(2)\n",
    "alphas_eu[0] = alphas_us[0]\n",
    "opt_eu = minimize(criterion_eu,[2.0,0.8],args=(sigma_c,phi,alphas_eu,yy_eu,mu_eu,[s_eu,h_eu]))\n",
    "alphas_eu[1] = opt_eu.x[0]\n",
    "p_eu = opt_eu.x[1]\n",
    "alphas_eu[1],p_eu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea21023f",
   "metadata": {},
   "source": [
    "We can check that we fit the targeted moments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "cb7573da",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0954833981139833, 0.9309709012192005)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ms_eu = [optm(sigma_c,phi,alphas_eu,i,p_eu,mu_eu,s_eu) for i in yy_eu]\n",
    "p_eu*np.mean(ms_eu)/np.mean(yy_eu), np.mean([prod(alphas_eu,m) for m in ms_eu])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b72224d5",
   "metadata": {},
   "source": [
    "The gradient is lower than in the U.S. but slighly larger than the target (1.06)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "7da1f53a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1453430675918532"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prod(alphas_eu,ms_eu[3])/prod(alphas_eu,ms_eu[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "271d4ee5",
   "metadata": {},
   "source": [
    "## Demand: Rand Experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "c9288990",
   "metadata": {},
   "outputs": [],
   "source": [
    "def focr(m,sigma_c,phi,alphas,y,p,mu,mom,mu0):\n",
    "\treturn mu*p*( (1-(1-mu0)*mom)*y - mu*p*m )**(-sigma_c) - alphas[0]*phi*np.exp(-alphas[1]-alphas[0]*m)\n",
    "\n",
    "def ur(m,sigma_c,phi,alphas,y,p,mu,mom,mu0):\n",
    "\treturn ( ((1-(1-mu0)*mom)*y - mu*p*m)**(1-sigma_c) )/(1-sigma_c) + phi*(1-np.exp(-alphas[1] - alphas[0]*m))\n",
    "\n",
    "def optmr(sigma_c,phi,alphas,y,p,mu,mom,mu0):\n",
    "\ttry:\n",
    "\t\tm = bisect(focr,0.0,y*0.99,args=(sigma_c,phi,alphas,y,p,mu,mom,mu0))\n",
    "\texcept:\n",
    "\t\tup  = ur(0.99*y,sigma_c,phi,alphas,y,p,mu,mom,mu0)\n",
    "\t\tlow = ur(0.0,sigma_c,phi,alphas,y,p,mu,mom,mu0)\n",
    "\t\tif up > low:\n",
    "\t\t\tm = 0.99*y\n",
    "\t\telse :\n",
    "\t\t\tm = 0.0\n",
    "\treturn m\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "b06244e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optmr(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us,mu_us) for i in yy_us]\n",
    "ms1 = [optmr(sigma_c,phi,alphas_us,i,p_us,mu_us*1.01,s_us,mu_us) for i in yy_us]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "576c3f59",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.6624542498869942"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (mu_us*1.01-mu_us)/mu_us\n",
    "\n",
    "((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1b2eb0f",
   "metadata": {},
   "source": [
    "## Demand : Price Elasticity in the US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "d9341d93",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optm(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us) for i in yy_us]\n",
    "ms1 = [optm(sigma_c,phi,alphas_us,i,p_us*1.01,mu_us,s_us) for i in yy_us]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "b28236ee",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.6624542498869936"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (p_us*1.01-p_us)/p_us\n",
    "\n",
    "((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "786586f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optm(sigma_c,phi,alphas_us,i,p_us*2,mu_us,s_us) for i in yy_us]\n",
    "ms1 = [optm(sigma_c,phi,alphas_us,i,p_us*2*1.01,mu_us,s_us) for i in yy_us]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "e6842fc2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.777901833570724"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (p_us*2*1.01-p_us*2)/(p_us*2)\n",
    "\n",
    "((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ee1b5c3",
   "metadata": {},
   "source": [
    "## Demand : Income Elasticity in the US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "4198d7fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optm(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us) for i in yy_us]\n",
    "ms1 = [optm(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us) for i in yy_us*1.01]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "892c5f2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2939210073015692"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (np.mean(1.01*yy_us)-np.mean(yy_us))/np.mean(yy_us)\n",
    "\n",
    "((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "9f47ac85",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optm(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us) for i in yy_us*2]\n",
    "ms1 = [optm(sigma_c,phi,alphas_us,i,p_us,mu_us,s_us) for i in yy_us*2*1.01]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "c4dc5d9b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8386611338909923"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (np.mean(1.01*yy_us*2)-np.mean(yy_us*2))/np.mean(yy_us*2)\n",
    "\n",
    "((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "183a75b8",
   "metadata": {},
   "source": [
    "## Demand : Price Elasticity in Europe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "f10b4443",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optm(sigma_c,phi,alphas_eu,i,p_eu,mu_eu,s_eu) for i in yy_eu]\n",
    "ms1 = [optm(sigma_c,phi,alphas_eu,i,p_eu*1.01,mu_eu,s_eu) for i in yy_eu]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "2d2c29ee",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.8166327822715324"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (p_eu*1.01-p_eu)/p_eu\n",
    "\n",
    "((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2469bccd",
   "metadata": {},
   "source": [
    "## Demand : Income Elasticity in Europe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "02f0e64d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ms0 = [optm(sigma_c,phi,alphas_eu,i,p_eu,mu_eu,s_eu) for i in yy_eu]\n",
    "ms1 = [optm(sigma_c,phi,alphas_eu,i,p_eu,mu_eu,s_eu) for i in yy_eu*1.01]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "744a8c18",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6104352012793126"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "level0 = np.mean(ms0)\n",
    "level1 = np.mean(ms1)\n",
    "dxox   = (np.mean(1.01*yy_eu)-np.mean(yy_eu))/np.mean(yy_eu)\n",
    "\n",
    "((level1 - level0)/level0)/dxox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3762e4c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "cf2a50979671a58939829e6829efb726aa5da11149213b77bd50351f899d04fb"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
